Chapter 6 Diversity analysis
6.1 Alpha diversity
# Calculate Hill numbers
richness <- genome_counts_filt %>%
column_to_rownames(var = "genome") %>%
dplyr::select(where(~ !all(. == 0))) %>%
hilldiv(., q = 0) %>%
t() %>%
as.data.frame() %>%
dplyr::rename(richness = 1) %>%
rownames_to_column(var = "sample")
neutral <- genome_counts_filt %>%
column_to_rownames(var = "genome") %>%
dplyr::select(where(~ !all(. == 0))) %>%
hilldiv(., q = 1) %>%
t() %>%
as.data.frame() %>%
dplyr::rename(neutral = 1) %>%
rownames_to_column(var = "sample")
phylogenetic <- genome_counts_filt %>%
column_to_rownames(var = "genome") %>%
dplyr::select(where(~ !all(. == 0))) %>%
hilldiv(., q = 1, tree = genome_tree) %>%
t() %>%
as.data.frame() %>%
dplyr::rename(phylogenetic = 1) %>%
rownames_to_column(var = "sample")
# Aggregate basal GIFT into elements
dist <- genome_gifts %>%
to.elements(., GIFT_db) %>%
traits2dist(., method = "gower")
functional <- genome_counts_filt %>%
column_to_rownames(var = "genome") %>%
dplyr::select(where(~ !all(. == 0))) %>%
hilldiv(., q = 1, dist = dist) %>%
t() %>%
as.data.frame() %>%
dplyr::rename(functional = 1) %>%
rownames_to_column(var = "sample") %>%
mutate(functional = if_else(is.nan(functional), 1, functional))
# Merge all metrics
alpha_div <- richness %>%
full_join(neutral, by = join_by(sample == sample)) %>%
full_join(phylogenetic, by = join_by(sample == sample)) %>%
full_join(functional, by = join_by(sample == sample))6.1.1 Wild samples
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="0_Wild") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
ggplot(aes(y = value, x = Population, group=Population, color=Population, fill=Population)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha=0.5) +
scale_color_manual(name="Population",
breaks=c("Cold_wet","Hot_dry"),
labels=c("Cold","Hot"),
values=c('#008080', "#d57d2c")) +
scale_fill_manual(name="Population",
breaks=c("Cold_wet","Hot_dry"),
labels=c("Cold","Hot"),
values=c('#00808050', "#d57d2c50")) +
facet_wrap(. ~ metric,scales = "free") +
coord_cartesian(xlim = c(1, NA)) +
stat_compare_means(size=3, label.x=.58) +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_text(size=10),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
) +
ylab("Alpha diversity")6.1.2 Acclimation samples
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="1_Acclimation") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
ggplot(aes(y = value, x = Population, group=Population, color=Population, fill=Population)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha=0.5) +
scale_color_manual(name="Population",
breaks=c("Cold_wet","Hot_dry"),
labels=c("Cold","Hot"),
values=c('#008080', "#d57d2c")) +
scale_fill_manual(name="Population",
breaks=c("Cold_wet","Hot_dry"),
labels=c("Cold","Hot"),
values=c('#00808050', "#d57d2c50")) +
facet_wrap(. ~ metric,scales = "free") +
coord_cartesian(xlim = c(1, NA)) +
stat_compare_means(size=3, label.x=.58) +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_text(size=10),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
) +
ylab("Alpha diversity")6.1.3 Antibiotics samples
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="2_Antibiotics") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
ggplot(aes(y = value, x = Population, group=Population, color=Population, fill=Population)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha=0.5) +
scale_color_manual(name="Population",
breaks=c("Cold_wet","Hot_dry"),
labels=c("Cold","Hot"),
values=c('#008080', "#d57d2c")) +
scale_fill_manual(name="Population",
breaks=c("Cold_wet","Hot_dry"),
labels=c("Cold","Hot"),
values=c('#00808050', "#d57d2c50")) +
facet_wrap(. ~ metric,scales = "free") +
coord_cartesian(xlim = c(1, NA)) +
stat_compare_means(size=3, label.x=.58) +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_text(size=10),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
) +
ylab("Alpha diversity")6.1.4 Transplant_1 samples
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="3_Transplant1") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha=0.5) +
scale_color_manual(name="Type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
values=c("#4477AA","#d57d2c", "#76b183")) +
scale_fill_manual(name="Type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
values=c("#4477AA50","#d57d2c50","#76b18350")) +
facet_wrap(. ~ metric,scales = "free") +
coord_cartesian(xlim = c(1, NA)) +
stat_compare_means(size=3, label.x=.7) +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_text(size=10),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
) +
ylab("Alpha diversity")6.1.5 Transplant_2 samples
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="4_Transplant2") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha=0.5) +
scale_color_manual(name="Type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
values=c("#4477AA","#d57d2c", "#76b183")) +
scale_fill_manual(name="Type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
values=c("#4477AA50","#d57d2c50","#76b18350")) +
facet_wrap(. ~ metric,scales = "free") +
coord_cartesian(xlim = c(1, NA)) +
stat_compare_means(size=3, label.x=.7) +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_text(size=10),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
) +
ylab("Alpha diversity")6.1.6 Post-Transplant_1 samples
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="5_Post-FMT1") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha=0.5) +
scale_color_manual(name="Type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
values=c("#4477AA","#d57d2c", "#76b183")) +
scale_fill_manual(name="Type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
values=c("#4477AA50","#d57d2c50","#76b18350")) +
facet_wrap(. ~ metric,scales = "free") +
coord_cartesian(xlim = c(1, NA)) +
stat_compare_means(size=3, label.x=.7) +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_text(size=10),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
) +
ylab("Alpha diversity")6.1.7 Post-Transplant_2 samples
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="6_Post-FMT2") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha=0.5) +
scale_color_manual(name="Type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
values=c("#4477AA","#d57d2c", "#76b183")) +
scale_fill_manual(name="Type",
breaks=c("Control","Hot_control", "Treatment"),
labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
values=c("#4477AA50","#d57d2c50","#76b18350")) +
facet_wrap(. ~ metric,scales = "free") +
coord_cartesian(xlim = c(1, NA)) +
stat_compare_means(size=3, label.x=.7) +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_text(size=10),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
) +
ylab("Alpha diversity")6.2 Beta diversity
beta_q0n <- genome_counts_filt %>%
column_to_rownames(., "genome") %>%
hillpair(., q = 0)
beta_q1n <- genome_counts_filt %>%
column_to_rownames(., "genome") %>%
hillpair(., q = 1)
beta_q1p <- genome_counts_filt %>%
column_to_rownames(., "genome") %>%
hillpair(., q = 1, tree = genome_tree)
beta_q1f <- genome_counts_filt %>%
column_to_rownames(., "genome") %>%
hillpair(., q = 1, dist = dist)6.3 Permanovas
6.3.1 1. Are the wild populations similar?
6.3.1.1 Wild: P.muralis vs P.liolepis
wild <- meta %>%
filter(time_point == "0_Wild")
# Create a temporary modified version of genome_counts_filt
temp_genome_counts <- transform(genome_counts_filt, row.names = genome_counts_filt$genome)
temp_genome_counts$genome <- NULL
wild.counts <- temp_genome_counts[, which(colnames(temp_genome_counts) %in% rownames(wild))]
identical(sort(colnames(wild.counts)), sort(as.character(rownames(wild))))
wild_nmds <- sample_metadata %>%
filter(time_point == "0_Wild")6.3.1.3 Richness
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.000012 0.000012 0.0012 999 0.979
Residuals 25 0.257281 0.010291
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Podarcis_liolepis Podarcis_muralis
Podarcis_liolepis 0.977
Podarcis_muralis 0.97302
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | 1.542719 | 0.2095041 | 6.625717 | 0.001 |
| Residual | 25 | 5.820951 | 0.7904959 | NA | NA |
| Total | 26 | 7.363669 | 1.0000000 | NA | NA |
6.3.1.4 Neutral
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.000048 0.0000476 0.0044 999 0.937
Residuals 25 0.270114 0.0108046
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Podarcis_liolepis Podarcis_muralis
Podarcis_liolepis 0.935
Podarcis_muralis 0.94763
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | 1.918266 | 0.2608511 | 8.822682 | 0.001 |
| Residual | 25 | 5.435610 | 0.7391489 | NA | NA |
| Total | 26 | 7.353876 | 1.0000000 | NA | NA |
6.3.1.5 Phylogenetic
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.03585 0.035847 2.4912 999 0.125
Residuals 25 0.35973 0.014389
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Podarcis_liolepis Podarcis_muralis
Podarcis_liolepis 0.125
Podarcis_muralis 0.12705
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | 0.3218613 | 0.2162815 | 6.899207 | 0.001 |
| Residual | 25 | 1.1662981 | 0.7837185 | NA | NA |
| Total | 26 | 1.4881594 | 1.0000000 | NA | NA |
6.3.1.6 Functional
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.018367 0.018367 1.5597 999 0.246
Residuals 25 0.294402 0.011776
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Podarcis_liolepis Podarcis_muralis
Podarcis_liolepis 0.234
Podarcis_muralis 0.22328
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | 0.0858578 | 0.172879 | 5.225323 | 0.051 |
| Residual | 25 | 0.4107775 | 0.827121 | NA | NA |
| Total | 26 | 0.4966352 | 1.000000 | NA | NA |
beta_q0n_nmds_wild <- beta_div_richness_wild$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace=FALSE) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(wild_nmds, by = join_by(sample == Tube_code))
beta_q1n_nmds_wild <- beta_div_neutral_wild$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace=FALSE) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(wild_nmds, by = join_by(sample == Tube_code))
beta_q1p_nmds_wild <- beta_div_phylo_wild$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(wild_nmds, by = join_by(sample == Tube_code))
beta_q1f_nmds_wild <- beta_div_func_wild$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(wild_nmds, by = join_by(sample == Tube_code))6.3.2 2. Effect of acclimation
accli <- meta %>%
filter(time_point == "1_Acclimation")
# Create a temporary modified version of genome_counts_filt
temp_genome_counts <- transform(genome_counts_filt, row.names = genome_counts_filt$genome)
temp_genome_counts$genome <- NULL
accli.counts <- temp_genome_counts[, which(colnames(temp_genome_counts) %in% rownames(accli))]
identical(sort(colnames(accli.counts)), sort(as.character(rownames(accli))))
accli_nmds <- sample_metadata %>%
filter(time_point == "1_Acclimation")6.3.2.2 Richness
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.11796 0.117959 12.963 999 0.003 **
Residuals 25 0.22748 0.009099
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Cold_wet Hot_dry
Cold_wet 0.003
Hot_dry 0.0013711
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| Population | 1 | 1.639807 | 0.179834 | 5.481634 | 0.001 |
| Residual | 25 | 7.478640 | 0.820166 | NA | NA |
| Total | 26 | 9.118447 | 1.000000 | NA | NA |
6.3.2.3 Neutral
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.07844 0.078443 5.2384 999 0.027 *
Residuals 25 0.37437 0.014975
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Cold_wet Hot_dry
Cold_wet 0.025
Hot_dry 0.030815
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| Population | 1 | 1.947003 | 0.2306127 | 7.493387 | 0.001 |
| Residual | 25 | 6.495736 | 0.7693873 | NA | NA |
| Total | 26 | 8.442739 | 1.0000000 | NA | NA |
6.3.2.4 Phylogenetic
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.06739 0.067395 2.9532 999 0.086 .
Residuals 25 0.57052 0.022821
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Cold_wet Hot_dry
Cold_wet 0.093
Hot_dry 0.098068
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| Population | 1 | 0.2441653 | 0.1224638 | 3.488854 | 0.019 |
| Residual | 25 | 1.7496100 | 0.8775362 | NA | NA |
| Total | 26 | 1.9937754 | 1.0000000 | NA | NA |
6.3.2.5 Functional
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.02496 0.024955 0.6729 999 0.41
Residuals 25 0.92714 0.037085
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Cold_wet Hot_dry
Cold_wet 0.39
Hot_dry 0.41979
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| Population | 1 | 0.0279454 | 0.0248037 | 0.6358634 | 0.447 |
| Residual | 25 | 1.0987171 | 0.9751963 | NA | NA |
| Total | 26 | 1.1266624 | 1.0000000 | NA | NA |
beta_q0n_nmds_accli <- beta_div_richness_accli$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace=FALSE) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(accli_nmds, by = join_by(sample == Tube_code))
beta_q1n_nmds_accli <- beta_div_neutral_accli$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace=FALSE) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(accli_nmds, by = join_by(sample == Tube_code))
beta_q1p_nmds_accli <- beta_div_phylo_accli$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(accli_nmds, by = join_by(sample == Tube_code))
beta_q1f_nmds_accli <- beta_div_func_accli$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(accli_nmds, by = join_by(sample == Tube_code))6.3.3 3. Comparison between Wild and Acclimation
accli1 <- meta %>%
filter(time_point == "0_Wild" | time_point == "1_Acclimation")
temp_genome_counts <- transform(genome_counts_filt, row.names = genome_counts_filt$genome)
temp_genome_counts$genome <- NULL
accli1.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(accli1))]
identical(sort(colnames(accli1.counts)),sort(as.character(rownames(accli1))))
accli1_nmds <- sample_metadata %>%
filter(time_point == "0_Wild" | time_point == "1_Acclimation")6.3.3.1 Number of samples used
[1] 54
beta_div_richness_accli1<-hillpair(data=accli1.counts, q=0)
beta_div_neutral_accli1<-hillpair(data=accli1.counts, q=1)
beta_div_phylo_accli1<-hillpair(data=accli1.counts, q=1, tree=genome_tree)
beta_div_func_accli1<-hillpair(data=accli1.counts, q=1, dist=dist)6.3.3.1.1 Richness
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.05014 0.050145 6.2252 999 0.014 *
Residuals 52 0.41886 0.008055
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
0_Wild 1_Acclimation
0_Wild 0.012
1_Acclimation 0.015808
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| time_point | 1 | 0.6172653 | 0.0360987 | 3.933397 | 0.001 |
| species | 1 | 2.8279677 | 0.1653842 | 18.020647 | 0.001 |
| individual | 25 | 9.5739861 | 0.5599025 | 2.440331 | 0.001 |
| Residual | 26 | 4.0801621 | 0.2386146 | NA | NA |
| Total | 53 | 17.0993812 | 1.0000000 | NA | NA |
6.3.3.1.2 Neutral
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.0199 0.0199035 2.1213 999 0.157
Residuals 52 0.4879 0.0093827
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
0_Wild 1_Acclimation
0_Wild 0.155
1_Acclimation 0.15128
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| time_point | 1 | 0.9050519 | 0.0541893 | 6.651487 | 0.001 |
| species | 1 | 3.3236300 | 0.1989999 | 24.426315 | 0.001 |
| individual | 25 | 8.9352276 | 0.5349902 | 2.626702 | 0.001 |
| Residual | 26 | 3.5377576 | 0.2118206 | NA | NA |
| Total | 53 | 16.7016671 | 1.0000000 | NA | NA |
6.3.3.1.3 Phylogenetic
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.01334 0.013340 0.6524 999 0.448
Residuals 52 1.06332 0.020449
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
0_Wild 1_Acclimation
0_Wild 0.453
1_Acclimation 0.42294
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| time_point | 1 | 0.2890434 | 0.0766494 | 7.532050 | 0.001 |
| species | 1 | 0.3508889 | 0.0930498 | 9.143655 | 0.001 |
| individual | 25 | 2.1332925 | 0.5657133 | 2.223620 | 0.004 |
| Residual | 26 | 0.9977533 | 0.2645874 | NA | NA |
| Total | 53 | 3.7709782 | 1.0000000 | NA | NA |
6.3.3.1.4 Functional
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.0123 0.012300 0.4817 999 0.53
Residuals 52 1.3277 0.025533
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
0_Wild 1_Acclimation
0_Wild 0.525
1_Acclimation 0.49073
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| time_point | 1 | 0.0448774 | 0.0269021 | 2.355512 | 0.172 |
| species | 1 | 0.0973005 | 0.0583275 | 5.107077 | 0.046 |
| individual | 25 | 1.0306426 | 0.6178264 | 2.163841 | 0.071 |
| Residual | 26 | 0.4953546 | 0.2969440 | NA | NA |
| Total | 53 | 1.6681751 | 1.0000000 | NA | NA |
beta_richness_nmds_accli1 <- beta_div_richness_accli1$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(accli1_nmds, by = c("sample" = "Tube_code"))
beta_neutral_nmds_accli1 <- beta_div_neutral_accli1$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(accli1_nmds, by = c("sample" = "Tube_code"))
beta_phylo_nmds_accli1 <- beta_div_phylo_accli1$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(accli1_nmds, by = join_by(sample == Tube_code))
beta_func_nmds_accli1 <- beta_div_func_accli1$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(accli1_nmds, by = join_by(sample == Tube_code))6.3.4 4. Do the antibiotics work?
6.3.4.1 Acclimation vs antibiotics
treat <- meta %>%
filter(time_point == "1_Acclimation" | time_point == "2_Antibiotics")
temp_genome_counts <- transform(genome_counts_filt, row.names = genome_counts_filt$genome)
temp_genome_counts$genome <- NULL
treat.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(treat))]
identical(sort(colnames(treat.counts)),sort(as.character(rownames(treat))))
treat_nmds <- sample_metadata %>%
filter(time_point == "1_Acclimation" | time_point == "2_Antibiotics")6.3.4.2 Number of samples used
[1] 50
beta_div_richness_treat<-hillpair(data=treat.counts, q=0)
beta_div_neutral_treat<-hillpair(data=treat.counts, q=1)
beta_div_phylo_treat<-hillpair(data=treat.counts, q=1, tree=genome_tree)
beta_div_func_treat<-hillpair(data=treat.counts, q=1, dist=dist)6.3.4.2.1 Richness
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.025318 0.0253178 6.021 999 0.014 *
Residuals 48 0.201837 0.0042049
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
1_Acclimation 2_Antibiotics
1_Acclimation 0.018
2_Antibiotics 0.017817
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| time_point | 1 | 1.888584 | 0.0949462 | 6.361098 | 0.001 |
| species | 1 | 2.117109 | 0.1064350 | 7.130814 | 0.001 |
| individual | 25 | 9.353701 | 0.4702455 | 1.260199 | 0.003 |
| Residual | 22 | 6.531709 | 0.3283734 | NA | NA |
| Total | 49 | 19.891103 | 1.0000000 | NA | NA |
6.3.4.2.2 Neutral
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.039587 0.039587 6.8387 999 0.009 **
Residuals 48 0.277854 0.005789
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
1_Acclimation 2_Antibiotics
1_Acclimation 0.01
2_Antibiotics 0.011886
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| time_point | 1 | 2.024181 | 0.1063620 | 9.051981 | 0.001 |
| species | 1 | 2.853103 | 0.1499183 | 12.758858 | 0.001 |
| individual | 25 | 9.234189 | 0.4852168 | 1.651783 | 0.001 |
| Residual | 22 | 4.919584 | 0.2585029 | NA | NA |
| Total | 49 | 19.031057 | 1.0000000 | NA | NA |
6.3.4.2.3 Phylogenetic
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.58372 0.58372 35.413 999 0.001 ***
Residuals 48 0.79119 0.01648
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
1_Acclimation 2_Antibiotics
1_Acclimation 0.001
2_Antibiotics 2.9795e-07
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| time_point | 1 | 1.8065206 | 0.2113909 | 18.636551 | 0.001 |
| species | 1 | 0.7903334 | 0.0924813 | 8.153292 | 0.001 |
| individual | 25 | 3.8164689 | 0.4465860 | 1.574869 | 0.007 |
| Residual | 22 | 2.1325541 | 0.2495419 | NA | NA |
| Total | 49 | 8.5458771 | 1.0000000 | NA | NA |
6.3.4.2.4 Functional
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.18591 0.185914 5.0679 999 0.027 *
Residuals 48 1.76088 0.036685
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
1_Acclimation 2_Antibiotics
1_Acclimation 0.03
2_Antibiotics 0.028989
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| time_point | 1 | 1.8020952 | 0.3750193 | 33.6195614 | 0.001 |
| species | 1 | 0.0031247 | 0.0006503 | 0.0582938 | 0.831 |
| individual | 25 | 1.8208629 | 0.3789249 | 1.3587875 | 0.192 |
| Residual | 22 | 1.1792568 | 0.2454055 | NA | NA |
| Total | 49 | 4.8053396 | 1.0000000 | NA | NA |
beta_richness_nmds_treat <- beta_div_richness_treat$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(treat_nmds, by = c("sample" = "Tube_code"))
beta_neutral_nmds_treat <- beta_div_neutral_treat$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(treat_nmds, by = c("sample" = "Tube_code"))
beta_phylo_nmds_treat <- beta_div_phylo_treat$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(treat_nmds, by = join_by(sample == Tube_code))
beta_func_nmds_treat <- beta_div_func_treat$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(treat_nmds, by = join_by(sample == Tube_code))6.3.5 5. Does the FMT work?
6.3.5.1 Comparison between FMT2 vs Post-FMT2
#Create newID to identify duplicated samples
transplants_metadata<-sample_metadata%>%
mutate(Tube_code=str_remove_all(Tube_code, "_a"))
transplants_metadata$newID <- paste(transplants_metadata$Tube_code, "_", transplants_metadata$individual)
transplant3<-transplants_metadata%>%
filter(time_point == "4_Transplant2" | time_point == "6_Post-FMT2")%>%
column_to_rownames("newID")
transplant3_nmds <- transplants_metadata %>%
filter(time_point == "4_Transplant2" | time_point == "6_Post-FMT2")
full_counts<-temp_genome_counts %>%
t()%>%
as.data.frame()%>%
rownames_to_column("Tube_code")%>%
full_join(transplants_metadata,by = join_by(Tube_code == Tube_code))
transplant3_counts<-full_counts %>%
filter(time_point == "4_Transplant2" | time_point == "6_Post-FMT2") %>%
subset(select=-c(315:324)) %>%
column_to_rownames("newID")%>%
subset(select=-c(1))%>%
t() %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric)
identical(sort(colnames(transplant3_counts)),sort(as.character(rownames(transplant3))))6.3.5.2 Number of samples used
[1] 49
beta_div_richness_transplant3<-hillpair(data=transplant3_counts, q=0)
beta_div_neutral_transplant3<-hillpair(data=transplant3_counts, q=1)
beta_div_phylo_transplant3<-hillpair(data=transplant3_counts, q=1, tree=genome_tree)
beta_div_func_transplant3<-hillpair(data=transplant3_counts, q=1, dist=dist)
#Arrange of metadata dataframe
transplant3_arrange<-transplant3[labels(beta_div_neutral_transplant3$S),]6.3.5.2.1 Richness
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | 1.180473 | 0.0855095 | 6.984555 | 0.001 |
| time_point | 1 | 0.860906 | 0.0623612 | 5.093759 | 0.001 |
| type | 1 | 1.459433 | 0.1057165 | 8.635089 | 0.001 |
| individual | 24 | 6.755100 | 0.4893170 | 1.665341 | 0.001 |
| Residual | 21 | 3.549250 | 0.2570959 | NA | NA |
| Total | 48 | 13.805162 | 1.0000000 | NA | NA |
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
1 Control vs Treatment 1 1.4169018 5.739828 0.15622903 0.001 0.003 *
2 Control vs Hot_control 1 2.0940966 8.509112 0.21005427 0.001 0.003 *
3 Treatment vs Hot_control 1 0.3004618 1.265034 0.04179854 0.144 0.432
6.3.5.2.2 Neutral
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | 1.2800927 | 0.0939787 | 8.796453 | 0.001 |
| time_point | 1 | 0.9350566 | 0.0686477 | 6.425458 | 0.001 |
| type | 1 | 1.9135997 | 0.1404879 | 13.149743 | 0.001 |
| individual | 24 | 6.4363516 | 0.4725281 | 1.842870 | 0.001 |
| Residual | 21 | 3.0559984 | 0.2243577 | NA | NA |
| Total | 48 | 13.6210990 | 1.0000000 | NA | NA |
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
1 Control vs Treatment 1 1.8758788 8.282671 0.21084796 0.001 0.003 *
2 Control vs Hot_control 1 2.4396317 10.635546 0.24945256 0.001 0.003 *
3 Treatment vs Hot_control 1 0.3158428 1.394345 0.04587515 0.115 0.345
6.3.5.2.3 Phylogenetic
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | 0.1400466 | 0.0952654 | 6.956436 | 0.001 |
| time_point | 1 | 0.1138047 | 0.0774145 | 5.652939 | 0.001 |
| type | 1 | 0.1432667 | 0.0974558 | 7.116383 | 0.001 |
| individual | 24 | 0.6501795 | 0.4422784 | 1.345663 | 0.046 |
| Residual | 21 | 0.4227709 | 0.2875859 | NA | NA |
| Total | 48 | 1.4700683 | 1.0000000 | NA | NA |
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
1 Control vs Treatment 1 0.14387705 5.735321 0.15612552 0.001 0.003 *
2 Control vs Hot_control 1 0.22715701 9.044894 0.22036587 0.001 0.003 *
3 Treatment vs Hot_control 1 0.04648319 1.704277 0.05550617 0.123 0.369
6.3.5.2.4 Functional
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | 0.0092808 | 0.0077189 | 0.4182529 | 0.528 |
| time_point | 1 | -0.0061674 | -0.0051295 | -0.2779456 | 0.883 |
| type | 1 | 0.0831052 | 0.0691191 | 3.7452726 | 0.101 |
| individual | 24 | 0.6501528 | 0.5407359 | 1.2208414 | 0.349 |
| Residual | 21 | 0.4659767 | 0.3875556 | NA | NA |
| Total | 48 | 1.2023481 | 1.0000000 | NA | NA |
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
1 Control vs Treatment 1 0.078539743 4.59293783 0.129040706 0.089 0.267
2 Control vs Hot_control 1 0.052468954 2.13675422 0.062593948 0.180 0.540
3 Treatment vs Hot_control 1 -0.002340352 -0.07432315 -0.002569452 0.857 1.000
beta_richness_nmds_transplant3 <- beta_div_richness_transplant3$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(transplant3_nmds, by = join_by(sample == newID))
beta_neutral_nmds_transplant3 <- beta_div_neutral_transplant3$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(transplant3_nmds, by = join_by(sample == newID))
beta_phylo_nmds_transplant3 <- beta_div_phylo_transplant3$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(transplant3_nmds, by = join_by(sample == newID))
beta_func_nmds_transplant3 <- beta_div_func_transplant3$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(transplant3_nmds, by = join_by(sample == newID))p0<-beta_richness_nmds_transplant3 %>%
group_by(individual) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
scale_color_manual(name="Type",
breaks=c("Control", "Hot_control", "Treatment"),
labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
values=c("#4477AA","#d57d2c","#76b183")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
theme_classic() +
theme(legend.position="none")
p1<-beta_neutral_nmds_transplant3 %>%
group_by(individual) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
scale_color_manual(name="Type",
breaks=c("Control", "Hot_control", "Treatment"),
labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
values=c("#4477AA","#d57d2c","#76b183")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
theme_classic() +
theme(legend.position="none")
p2<-beta_phylo_nmds_transplant3 %>%
group_by(individual) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
scale_color_manual(name="Type",
breaks=c("Control", "Hot_control", "Treatment"),
labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
values=c("#4477AA","#d57d2c","#76b183")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
theme_classic() +
theme(legend.position="none")
p3<-beta_func_nmds_transplant3 %>%
group_by(individual) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
scale_color_manual(name="Type",
breaks=c("Control", "Hot_control", "Treatment"),
labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
values=c("#4477AA","#d57d2c","#76b183")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y= element_blank (), x="NMDS1 \n Functional beta diversity") +
theme_classic()+
theme(legend.position="none")6.3.5.3 Comparison between the different experimental time points
6.3.5.3.1 Comparison againts Wild samples
transplants_metadata<-sample_metadata%>%
mutate(Tube_code=str_remove_all(Tube_code, "_a"))
transplants_metadata$newID <- paste(transplants_metadata$Tube_code, "_", transplants_metadata$individual)
transplant4<-transplants_metadata%>%
filter(time_point == "4_Transplant2" | time_point == "6_Post-FMT2" |time_point == "0_Wild")%>%
column_to_rownames("newID")
transplant4_pairwise <- transplants_metadata %>%
filter(time_point == "4_Transplant2" | time_point == "6_Post-FMT2"|time_point == "0_Wild")
full_counts<-temp_genome_counts %>%
t()%>%
as.data.frame()%>%
rownames_to_column("Tube_code")%>%
full_join(transplants_metadata,by = join_by(Tube_code == Tube_code))
transplant4_counts<-full_counts %>%
filter(time_point == "4_Transplant2" | time_point == "6_Post-FMT2"|time_point == "0_Wild") %>%
subset(select=-c(315:324)) %>%
column_to_rownames("newID")%>%
subset(select=-c(1))%>%
t() %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric)
identical(sort(colnames(transplant4_counts)),sort(as.character(rownames(transplant4))))
beta_div_neutral_transplant4<-hillpair(data=transplant4_counts, q=1)The estimated time for calculating the 2850 pairwise combinations is 13 seconds.
transplant4_beta_div_neutral<-as.matrix(beta_div_neutral_transplant4$S) #convert the pairwise comparisons into a matrix
transplant4_beta_div_neutral_df <- as.data.frame(as.table(as.matrix(transplant4_beta_div_neutral))) #convert the matrix to a data frame with pairwise valuestransplant4_beta_div_neutral_df <- transplant4_beta_div_neutral_df %>% #remove the duplicated pairwise comparisons
filter(Var1 != Var2) %>%
mutate(pair = pmap_chr(list(Var1, Var2), ~paste(sort(c(..1, ..2)), collapse = "_"))) %>%
distinct(pair, .keep_all = TRUE) %>%
select(Sample_1 = Var1, Sample_2 = Var2, Distance = Freq) %>%
arrange(Sample_1, Sample_2)
transplant4_beta_div_neutral_df <- transplant4_beta_div_neutral_df %>% #keep only the pairwise comparisons between the same individual
mutate(Sample_1_part = sub("^[^_]*_", "", Sample_1),
Sample_2_part = sub("^[^_]*_", "", Sample_2)) %>%
filter(Sample_1_part == Sample_2_part)
#%>%
#select(Sample_1, Sample_2, Distance)
transplant4_beta_div_neutral_met<-transplant4_beta_div_neutral_df %>% #merge with the metadata associated to the pairwise comparisons for Sample_1
inner_join(transplant4_pairwise, by = join_by(Sample_1 == "newID"))
# Extract the part before the first `_` for Sample_1 and Sample_2
transplant4_beta_div_neutral_met$Sample_1_part <- sub("_.*", "", transplant4_beta_div_neutral_met$Sample_1)
# Create a named vector for matching Tube_code with Time_point
time_point_map <- setNames(transplant4_beta_div_neutral_met$time_point, transplant4_beta_div_neutral_met$Tube_code)
# Function to replace the part before the first _ with the corresponding Time_point
replace_with_time_point <- function(sample_part, tube_code, time_point_map) {
if (tube_code %in% names(time_point_map)) {
return(time_point_map[tube_code])
} else {
return(sample_part)
}
}
# Apply the function to Sample_1_part
transplant4_beta_div_neutral_met$Sample_1_part <- mapply(replace_with_time_point,
transplant4_beta_div_neutral_met$Sample_1_part,
transplant4_beta_div_neutral_met$Tube_code,
list(time_point_map))
transplant4_beta_div_neutral_met<-transplant4_beta_div_neutral_met %>%
subset(select=-c(6:16))
transplant4_beta_div_neutral_met$Sample_2_part <- sub("_.*", "", transplant4_beta_div_neutral_met$Sample_2)
transplant4_beta_div_neutral_met<-transplant4_beta_div_neutral_met %>% #merge with the metadata associated to the pairwise comparisons for Sample_2
inner_join(transplant4_pairwise, by = join_by(Sample_2 == "newID"))
# Create a named vector again for matching Tube_code with Time_point
time_point_map <- setNames(transplant4_beta_div_neutral_met$time_point, transplant4_beta_div_neutral_met$Tube_code)
replace_with_time_point <- function(sample_part, tube_code, time_point_map) {
if (tube_code %in% names(time_point_map)) {
return(time_point_map[tube_code])
} else {
return(sample_part)
}
}
# Apply the function to Sample_2_part
transplant4_beta_div_neutral_met$Sample_2_part <- mapply(replace_with_time_point,
transplant4_beta_div_neutral_met$Sample_2_part,
transplant4_beta_div_neutral_met$Tube_code,
list(time_point_map))
# Create a new variable to plot the distances between the time_points
transplant4_beta_div_neutral_met$comparison <- paste(transplant4_beta_div_neutral_met$Sample_1_part, "-",transplant4_beta_div_neutral_met$Sample_2_part)##plot the differences between wild and post-transplant and transplant and post-transplant
transplant4_beta_div_neutral_met %>%
ggplot(aes(x=type, y=Distance, color=type, fill=type))+
geom_boxplot() +
scale_color_manual(name="Type",
breaks=c("Control", "Hot_control", "Treatment"),
labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
values=c("#4477AA","#d57d2c","#76b183")) +
scale_fill_manual(name="Type",
breaks=c("Control", "Hot_control", "Treatment"),
labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
values=c("#4477AA50","#d57d2c50","#76b18350")) +
scale_x_discrete(labels=c("Control" = "Cold-Cold", "Hot_control" = "Hot-Hot", "Treatment" = "Cold-Hot")) +
theme_minimal() +
facet_wrap(~comparison, labeller = labeller(comparison = c("0_Wild - 6_Post-FMT2" = "Wild-Treatment", "0_Wild - 4_Transplant2" = "Wild-Transplant", "6_Post-FMT2 - 4_Transplant2"="Treatment-Transplant"))) +
theme(
axis.text.x = element_text(angle = 45, hjust = 1))6.3.5.3.2 Comparison againts Acclimation samples
#Create newID to identify duplicated samples
transplants_metadata<-sample_metadata%>%
mutate(Tube_code=str_remove_all(Tube_code, "_a"))
transplants_metadata$newID <- paste(transplants_metadata$Tube_code, "_", transplants_metadata$individual)
transplant5<-transplants_metadata%>%
filter(time_point == "4_Transplant2" | time_point == "6_Post-FMT2" |time_point == "1_Acclimation")%>%
column_to_rownames("newID")
transplant5_pairwise <- transplants_metadata %>%
filter(time_point == "4_Transplant2" | time_point == "6_Post-FMT2"|time_point == "1_Acclimation")
full_counts_new<-temp_genome_counts %>%
t()%>%
as.data.frame()%>%
rownames_to_column("Tube_code")%>%
full_join(transplants_metadata,by = join_by(Tube_code == Tube_code))
transplant5_counts<-full_counts %>%
filter(time_point == "4_Transplant2" | time_point == "6_Post-FMT2"|time_point == "1_Acclimation") %>%
subset(select=-c(315:324)) %>%
column_to_rownames("newID")%>%
subset(select=-c(1))%>%
t() %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric)
identical(sort(colnames(transplant5_counts)),sort(as.character(rownames(transplant5))))
beta_div_neutral_transplant5<-hillpair(data=transplant5_counts, q=1)The estimated time for calculating the 2850 pairwise combinations is 13 seconds.
transplant5_beta_div_neutral<-as.matrix(beta_div_neutral_transplant5$S) #convert the pairwise comparisons into a matrix
transplant5_beta_div_neutral_df <- as.data.frame(as.table(as.matrix(transplant5_beta_div_neutral))) #convert the matrix to a data frame with pairwise valuestransplant5_beta_div_neutral_df <- transplant5_beta_div_neutral_df %>% #remove the duplicated pairwise comparisons
filter(Var1 != Var2) %>%
mutate(pair = pmap_chr(list(Var1, Var2), ~paste(sort(c(..1, ..2)), collapse = "_"))) %>%
distinct(pair, .keep_all = TRUE) %>%
select(Sample_1 = Var1, Sample_2 = Var2, Distance = Freq) %>%
arrange(Sample_1, Sample_2)
transplant5_beta_div_neutral_df <- transplant5_beta_div_neutral_df %>% #keep only the pairwise comparisons between the same individual
mutate(Sample_1_part = sub("^[^_]*_", "", Sample_1),
Sample_2_part = sub("^[^_]*_", "", Sample_2)) %>%
filter(Sample_1_part == Sample_2_part)
#%>%
#select(Sample_1, Sample_2, Distance)
transplant5_beta_div_neutral_met<-transplant5_beta_div_neutral_df %>% #merge with the metadata associated to the pairwise comparisons for Sample_1
inner_join(transplant5_pairwise, by = join_by(Sample_1 == "newID"))
# Extract the part before the first `_` for Sample_1 and Sample_2
transplant5_beta_div_neutral_met$Sample_1_part <- sub("_.*", "", transplant5_beta_div_neutral_met$Sample_1)
# Create a named vector for matching Tube_code with Time_point
time_point_map_new <- setNames(transplant5_beta_div_neutral_met$time_point, transplant5_beta_div_neutral_met$Tube_code)
# Function to replace the part before the first _ with the corresponding Time_point
replace_with_time_point_new <- function(sample_part, tube_code, time_point_map_new) {
if (tube_code %in% names(time_point_map_new)) {
return(time_point_map_new[tube_code])
} else {
return(sample_part)
}
}
# Apply the function to Sample_1_part
transplant5_beta_div_neutral_met$Sample_1_part <- mapply(replace_with_time_point_new,
transplant5_beta_div_neutral_met$Sample_1_part,
transplant5_beta_div_neutral_met$Tube_code,
list(time_point_map_new))
transplant5_beta_div_neutral_met<-transplant5_beta_div_neutral_met %>%
subset(select=-c(6:16))
transplant5_beta_div_neutral_met$Sample_2_part <- sub("_.*", "", transplant5_beta_div_neutral_met$Sample_2)
transplant5_beta_div_neutral_met<-transplant5_beta_div_neutral_met %>% #merge with the metadata associated to the pairwise comparisons for Sample_2
inner_join(transplant5_pairwise, by = join_by(Sample_2 == "newID"))
# Create a named vector again for matching Tube_code with Time_point
time_point_map_new <- setNames(transplant5_beta_div_neutral_met$time_point, transplant5_beta_div_neutral_met$Tube_code)
replace_with_time_point_new <- function(sample_part, tube_code, time_point_map_new) {
if (tube_code %in% names(time_point_map_new)) {
return(time_point_map_new[tube_code])
} else {
return(sample_part)
}
}
# Apply the function to Sample_2_part
transplant5_beta_div_neutral_met$Sample_2_part <- mapply(replace_with_time_point_new,
transplant5_beta_div_neutral_met$Sample_2_part,
transplant5_beta_div_neutral_met$Tube_code,
list(time_point_map_new))
# Create a new variable to plot the distances between the time_points
transplant5_beta_div_neutral_met$comparison <- paste(transplant5_beta_div_neutral_met$Sample_1_part, "-",transplant5_beta_div_neutral_met$Sample_2_part)
##plot the differences between acclimation and post-transplant and transplant and post-transplant
# Reorder the 'comparison' factor
transplant5_beta_div_neutral_met$comparison <- factor(
transplant5_beta_div_neutral_met$comparison,
levels = c("4_Transplant2 - 1_Acclimation", # Acclimation-Transplant first
"6_Post-FMT2 - 4_Transplant2", # Transplant-Treatment second
"6_Post-FMT2 - 1_Acclimation") # Treatment-Acclimation last
)# Create the plot with the specified facet order
transplant5_beta_div_neutral_met %>%
ggplot(aes(x=type, y=Distance, color=type, fill=type)) +
geom_boxplot() +
scale_color_manual(name="Type",
breaks=c("Control", "Hot_control", "Treatment"),
labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
values=c("#4477AA","#d57d2c","#76b183")) +
scale_fill_manual(name="Type",
breaks=c("Control", "Hot_control", "Treatment"),
labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
values=c("#4477AA50","#d57d2c50","#76b18350")) +
scale_x_discrete(labels=c("Control" = "Cold-Cold", "Hot_control" = "Hot-Hot", "Treatment" = "Cold-Hot")) +
theme_minimal() +
facet_wrap(~comparison,
labeller = labeller(comparison = c(
"4_Transplant2 - 1_Acclimation" = "Acclimation-Transplant",
"6_Post-FMT2 - 4_Transplant2" = "Transplant-Treatment",
"6_Post-FMT2 - 1_Acclimation" = "Treatment-Acclimation"
)),
nrow = 1, ncol = 3) + # Single row with 3 columns
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
ylab("Beta diversity dissimilarity distance") +
xlab("Type")6.3.6 7. Are there differences between the control and the treatment group?
6.3.6.2 Number of samples used
[1] 26
beta_div_richness_post1<-hillpair(data=post1.counts, q=0)
beta_div_neutral_post1<-hillpair(data=post1.counts, q=1)
beta_div_phylo_post1<-hillpair(data=post1.counts, q=1, tree=genome_tree)
beta_div_func_post1<-hillpair(data=post1.counts, q=1, dist=dist)
#Arrange of metadata dataframe
post1_arrange<-post1[labels(beta_div_neutral_post1$S),]6.3.6.2.1 Richness
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.017675 0.0088373 2.3825 999 0.095 .
Residuals 23 0.085312 0.0037092
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Control Hot_control Treatment
Control 0.0030000 0.647
Hot_control 0.0068795 0.215
Treatment 0.6248469 0.2084296
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | 0.6340254 | 0.0768024 | 2.065607 | 0.004 |
| type | 1 | 0.5615418 | 0.0680222 | 1.829461 | 0.009 |
| Residual | 23 | 7.0597099 | 0.8551754 | NA | NA |
| Total | 25 | 8.2552771 | 1.0000000 | NA | NA |
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
1 Control vs Treatment 1 0.5615418 1.729004 0.1033537 0.020 0.060
2 Control vs Hot_control 1 0.8438429 2.793772 0.1486541 0.001 0.003 *
3 Treatment vs Hot_control 1 0.3734921 1.268929 0.0779971 0.103 0.309
6.3.6.2.2 Neutral
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.011001 0.0055005 0.6303 999 0.54
Residuals 23 0.200714 0.0087267
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Control Hot_control Treatment
Control 0.21300 0.961
Hot_control 0.21166 0.448
Treatment 0.95468 0.43604
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | 0.7907904 | 0.1076445 | 3.056657 | 0.001 |
| type | 1 | 0.6051778 | 0.0823784 | 2.339205 | 0.010 |
| Residual | 23 | 5.9503501 | 0.8099772 | NA | NA |
| Total | 25 | 7.3463184 | 1.0000000 | NA | NA |
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
1 Control vs Treatment 1 0.6051778 2.250849 0.13047758 0.017 0.051
2 Control vs Hot_control 1 1.0528902 4.143637 0.20570451 0.001 0.003 *
3 Treatment vs Hot_control 1 0.4150076 1.637268 0.09840968 0.044 0.132
6.3.6.2.3 Phylogenetic
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.00440 0.0021994 0.1369 999 0.926
Residuals 23 0.36941 0.0160614
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Control Hot_control Treatment
Control 0.94500 0.719
Hot_control 0.91505 0.775
Treatment 0.63312 0.73046
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | 0.0560850 | 0.0531376 | 1.3149967 | 0.262 |
| type | 1 | 0.0184254 | 0.0174571 | 0.4320099 | 0.791 |
| Residual | 23 | 0.9809570 | 0.9294053 | NA | NA |
| Total | 25 | 1.0554673 | 1.0000000 | NA | NA |
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
1 Control vs Treatment 1 0.01842535 0.4144162 0.02688498 0.782 1.000
2 Control vs Hot_control 1 0.05987967 1.7387847 0.09802164 0.106 0.318
3 Treatment vs Hot_control 1 0.03212966 0.6477782 0.04139746 0.690 1.000
6.3.6.2.4 Functional
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.00400 0.0020014 0.145 999 0.873
Residuals 23 0.31753 0.0138057
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Control Hot_control Treatment
Control 0.61200 0.745
Hot_control 0.59817 0.849
Treatment 0.75141 0.83718
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| species | 1 | 0.0024979 | 0.0033024 | 0.0900845 | 0.625 |
| type | 1 | 0.1161466 | 0.1535542 | 4.1887855 | 0.076 |
| Residual | 23 | 0.6377435 | 0.8431434 | NA | NA |
| Total | 25 | 0.7563879 | 1.0000000 | NA | NA |
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
1 Control vs Treatment 1 0.11614656 4.724791 0.23953568 0.073 0.219
2 Control vs Hot_control 1 0.05000930 1.704826 0.09629160 0.229 0.687
3 Treatment vs Hot_control 1 0.01235859 0.423812 0.02747777 0.471 1.000
beta_richness_nmds_post1 <- beta_div_richness_post1$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(post1_nmds, by = join_by(sample == Tube_code))
beta_neutral_nmds_post1 <- beta_div_neutral_post1$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(post1_nmds, by = join_by(sample == Tube_code))
beta_phylogenetic_nmds_post1 <- beta_div_phylo_post1$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(post1_nmds, by = join_by(sample == Tube_code))
beta_functional_nmds_post1 <- beta_div_func_post1$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(post1_nmds, by = join_by(sample == Tube_code))p0<-beta_richness_nmds_post1 %>%
group_by(type) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
scale_color_manual(name="Type",
breaks=c("Control", "Hot_control", "Treatment"),
labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
values=c("#4477AA","#d57d2c","#76b183")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
theme_classic() +
theme(legend.position="none")
p1<-beta_neutral_nmds_post1 %>%
group_by(type) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
scale_color_manual(name="Type",
breaks=c("Control", "Hot_control", "Treatment"),
labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
values=c("#4477AA","#d57d2c","#76b183")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
theme_classic() +
theme(legend.position="none")
p2<-beta_phylogenetic_nmds_post1 %>%
group_by(type) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
scale_color_manual(name="Type",
breaks=c("Control", "Hot_control", "Treatment"),
labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
values=c("#4477AA","#d57d2c","#76b183")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
theme_classic() +
theme(legend.position="none")
p3<-beta_functional_nmds_post1 %>%
group_by(type) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
scale_color_manual(name="Type",
breaks=c("Control", "Hot_control", "Treatment"),
labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
values=c("#4477AA","#d57d2c","#76b183")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y= element_blank (), x="NMDS1 \n Functional beta diversity") +
theme_classic()+
theme(legend.position="none")6.3.6.4 Number of samples used
[1] 27
beta_div_richness_post2<-hillpair(data=post2.counts, q=0)
beta_div_neutral_post2<-hillpair(data=post2.counts, q=1)
beta_div_phylo_post2<-hillpair(data=post2.counts, q=1, tree=genome_tree)
beta_div_func_post2<-hillpair(data=post2.counts, q=1, dist=dist)
#Arrange of metadata dataframe
post2_arrange<-post2[labels(beta_div_neutral_post2$S),]6.3.6.4.1 Richness
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.002011 0.0010056 0.1982 999 0.818
Residuals 24 0.121775 0.0050740
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Control Hot_control Treatment
Control 0.71700 0.807
Hot_control 0.67789 0.571
Treatment 0.79246 0.59820
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| type | 2 | 1.504341 | 0.1967776 | 2.939822 | 0.001 |
| Residual | 24 | 6.140538 | 0.8032224 | NA | NA |
| Total | 26 | 7.644879 | 1.0000000 | NA | NA |
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
1 Treatment vs Control 1 0.6463814 2.560441 0.1379515 0.001 0.003 *
2 Treatment vs Hot_control 1 0.4796256 1.916520 0.1069694 0.003 0.009 *
3 Control vs Hot_control 1 1.1305044 4.268317 0.2105906 0.001 0.003 *
6.3.6.4.2 Neutral
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.008262 0.0041311 0.8024 999 0.479
Residuals 24 0.123559 0.0051483
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Control Hot_control Treatment
Control 0.42400 0.682
Hot_control 0.44675 0.256
Treatment 0.65989 0.25095
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| type | 2 | 1.923807 | 0.2603795 | 4.224537 | 0.001 |
| Residual | 24 | 5.464666 | 0.7396205 | NA | NA |
| Total | 26 | 7.388473 | 1.0000000 | NA | NA |
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
1 Treatment vs Control 1 1.0227481 4.648335 0.2251191 0.001 0.003 *
2 Treatment vs Hot_control 1 0.5010202 2.206532 0.1211945 0.002 0.006 *
3 Control vs Hot_control 1 1.3619424 5.771031 0.2650785 0.001 0.003 *
6.3.6.4.3 Phylogenetic
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.000407 0.0002034 0.0487 999 0.969
Residuals 24 0.100305 0.0041794
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Control Hot_control Treatment
Control 0.95600 0.857
Hot_control 0.93765 0.775
Treatment 0.83933 0.76015
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| type | 2 | 0.1594363 | 0.2042241 | 3.079623 | 0.001 |
| Residual | 24 | 0.6212564 | 0.7957759 | NA | NA |
| Total | 26 | 0.7806927 | 1.0000000 | NA | NA |
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
1 Treatment vs Control 1 0.05927454 2.382025 0.1295845 0.033 0.099
2 Treatment vs Hot_control 1 0.06906280 2.722460 0.1454115 0.004 0.012 .
3 Control vs Hot_control 1 0.11081709 4.043656 0.2017424 0.002 0.006 *
6.3.6.4.4 Functional
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.01126 0.0056302 0.2861 999 0.811
Residuals 24 0.47233 0.0196806
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Control Hot_control Treatment
Control 0.53100 0.680
Hot_control 0.48255 0.813
Treatment 0.60116 0.75643
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| type | 2 | -0.0038724 | -0.0056213 | -0.0670788 | 0.92 |
| Residual | 24 | 0.6927468 | 1.0056213 | NA | NA |
| Total | 26 | 0.6888744 | 1.0000000 | NA | NA |
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
1 Treatment vs Control 1 -0.008527330 -0.46290555 -0.029793572 0.856 1
2 Treatment vs Hot_control 1 -0.001648721 -0.04717131 -0.002956924 0.908 1
3 Control vs Hot_control 1 0.004367477 0.13147026 0.008149924 0.679 1
beta_richness_nmds_post2 <- beta_div_richness_post2$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(post2_nmds, by = join_by(sample == Tube_code))
beta_neutral_nmds_post2 <- beta_div_neutral_post2$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(post2_nmds, by = join_by(sample == Tube_code))
beta_phylogenetic_nmds_post2 <- beta_div_phylo_post2$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(post2_nmds, by = join_by(sample == Tube_code))
beta_functional_nmds_post2 <- beta_div_func_post2$S %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(post2_nmds, by = join_by(sample == Tube_code))p0<-beta_richness_nmds_post2 %>%
group_by(type) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
scale_color_manual(name="Type",
breaks=c("Control", "Hot_control", "Treatment"),
labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
values=c("#4477AA","#d57d2c","#76b183")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
theme_classic() +
theme(legend.position="none")
p1<-beta_neutral_nmds_post2 %>%
group_by(type) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
scale_color_manual(name="Type",
breaks=c("Control", "Hot_control", "Treatment"),
labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
values=c("#4477AA","#d57d2c","#76b183")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
theme_classic() +
theme(legend.position="none")
p2<-beta_phylogenetic_nmds_post2 %>%
group_by(type) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
scale_color_manual(name="Type",
breaks=c("Control", "Hot_control", "Treatment"),
labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
values=c("#4477AA","#d57d2c","#76b183")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
theme_classic() +
theme(legend.position="none")
p3<-beta_functional_nmds_post2 %>%
group_by(type) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
scale_color_manual(name="Type",
breaks=c("Control", "Hot_control", "Treatment"),
labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
values=c("#4477AA","#d57d2c","#76b183")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y= element_blank (), x="NMDS1 \n Functional beta diversity") +
theme_classic()+
theme(legend.position="none")